Mixed-Effects Models

Repeated Measures Model Dataset: Mohlman, Jessica L. et al. (2019), Data from: Nonconsumptive effects of hunting on a nontarget game bird, Dryad Dataset

Disturbance through human hunting activity can significantly impact prey species through both consumptive and nonconsumptive effects. The nonconsumptive effects of rabbit hunting on Northern Bobwhite (Colinus virginianus; hereafter, bobwhite) could cause an increased perceived risk of predation by bobwhite during rabbit hunting events may elicit anti-predator responses, such as reduced movement away from the safety of cover.

bobwhite <- read.csv('bobwhite3.csv')
bobwhite$ID <- as.factor(bobwhite$ID) #make ID a factor
p1<- ggplot(bobwhite, aes(x=HuntDay, y=HW_Dist, group=ID, color=ID, shape=ID)) + 
  geom_point(size=4, alpha=0.6, position = position_dodge2(width=.33, preserve = "total")) +
  scale_y_continuous() +
  #geom_line() +
  geom_smooth(method = "lm", se = FALSE) +
  labs(title="Risk Behavior in Bobwhite During Hunting Season", x= "Hunting Season Species", y = "Distance from Hardwood Forest Cover (meters)")+ 
  theme_bw()+
  scale_color_brewer(palette = "BrBG")
p1

## `geom_smooth()` using formula 'y ~ x'

# ID: Unique ID given to each bobwhite covey tracked in chronological order.
# HuntDay: Denotes if it was a “Rabbit” or “Quail” (bobwhite) scheduled hunt day. 
# HW_Dist: Distance in meters a bobwhite covey was from hardwood habitat.

bobwhite_means <- bobwhite %>%
  group_by(HuntDay) %>%
  summarise(mean_HW_Dist=mean(HW_Dist),
            se_HW_Dist=sd(HW_Dist)/sqrt(n()))
bobwhite_means

## # A tibble: 2 × 3
##   HuntDay mean_HW_Dist se_HW_Dist
##   <chr>          <dbl>      <dbl>
## 1 Quail           22.3       4.29
## 2 Rabbit          57.0       5.30

mixed_bobwhite <- lmer(HW_Dist~(HuntDay*ID)+(1|ID), data = bobwhite)

## Warning in as_lmerModLT(model, devfun): Model may not have converged with 1
## eigenvalue close to zero: 5.5e-09

anova(mixed_bobwhite) 

## Type III Analysis of Variance Table with Satterthwaite's method
##             Sum Sq Mean Sq NumDF DenDF F value   Pr(>F)   
## HuntDay    19204.9 19204.9     1   188 11.0107 0.001088 **
## ID          3168.5   633.7     5   188  0.3633 0.873151   
## HuntDay:ID 23341.2  4668.2     5   188  2.6764 0.023103 * 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

performance::check_model(mixed_bobwhite)

bob_means <- emmeans(mixed_bobwhite, "HuntDay")

## Cannot use mode = "kenward-roger" because *pbkrtest* package is not installed

## NOTE: Results may be misleading due to involvement in interactions

bob_means

##  HuntDay emmean   SE  df lower.CL upper.CL
##  Quail     23.1 25.3 188   -26.92     73.1
##  Rabbit    45.1 25.1 188    -4.31     94.5
## 
## Results are averaged over the levels of: ID 
## Degrees-of-freedom method: satterthwaite 
## Confidence level used: 0.95

Because there is an unequal number of points per day for each individual covey, we can check the accuracy of this model by using emmeans:

em_bob <- as.data.frame(bob_means)
em_bob

##   HuntDay   emmean       SE  df   lower.CL upper.CL
## 1   Quail 23.08173 25.34776 188 -26.920856 73.08431
## 2  Rabbit 45.11574 25.05572 188  -4.310749 94.54224

bobwhite_means

## # A tibble: 2 × 3
##   HuntDay mean_HW_Dist se_HW_Dist
##   <chr>          <dbl>      <dbl>
## 1 Quail           22.3       4.29
## 2 Rabbit          57.0       5.30

em_bob$HuntDay <- factor(em_bob$HuntDay, levels= c("Quail","Rabbit"))

ggplot(em_bob, aes(x=HuntDay, y=emmean)) + 
  geom_point(size=5, color="#87b8d3") +
  geom_errorbar(aes(ymin=emmean-SE, ymax=emmean+SE), width=.2, color="#6699CC") +
  geom_point(size=5, data=bobwhite_means, x=bobwhite_means$HuntDay, y=bobwhite_means$mean_HW_Dist, color = "#336699") +
  theme(axis.text.x = ggtext::element_markdown(color = "tan4", size = 12)) +
  scale_x_discrete(labels = labels) +
  theme(plot.caption=element_text(size=9, hjust=0, margin=margin(15,0,0,0)))+
  theme_bw()+
  labs(title="Comparing raw means to emmeans", caption="Light blue = raw means, Dark blue = adjusted means", x= "Hunting Day (1=Quail, 2=Rabbit)", y = "Average distance from hardwood cover (m)") 

\(~\) \(~\)

Nested Hierarchical Model

Dataset: Snijders, Lysanne; van Oers, Kees; Naguib, Marc (2017), Data from: Sex-specific responses to territorial intrusions in a communication network: evidence from radio-tagged great tits, Dryad Dataset

Great Tits (Parus major)

Advertisement signaling is usually linked to intersexual selection and intrasexual competition and thus is a key component of a species’ ecology. Using a novel spatial tracking system, the authors tested whether or not the spatial behavior of male and female great tits (Parus major) changes in relation to the response of a territorial male neighbor to an intruder. They tracked the spatial behavior of male and female great tits (N = 13), 1 hr before and 1 hr after simulating territory intrusions. The individual bird is a random effect, and the sex of the bird and measurements taken before and after are fixed effects.

tit <- read.csv("great_tits.csv")
tit$sex <- as.factor(tit$sex)
ggplot(tit, aes(sex, dist_m, colour = as.factor(id), shape=as.factor(measure))) + 
  geom_jitter(width =0.15, size=5, alpha=0.6)+
  ylab ("Response Distance from Playback Location (meters)") +
  xlab ("Sex") +
  annotate("text", x = 2, y = 76, label = "13 birds") +
  annotate("text", x = 2, y = 72, label = "2 measurements per bird") +
  annotate("text", x = 2, y = 67, label = "26 total measurements", size=4, color="blue")+
  labs(title="Sex-specific Responses to Territorial Intrusions", caption="1 = female
       2 = male")+
  scale_shape_discrete(name= "Measurements 
  Before / After")+
  theme_bw()

In this study, 13 individuals were measured 2 times, 1 hour before playback and 1 hour after. The individuals are separated by sex to look at the different distances between females and males. Distance on the Y axis is measured in meters, with distance being how many meters away a bird was detected before and after a playback recording was played of a territorial male’s song. Based on the measurements of distance before and after playback, it looks like there isn’t really a pattern of whether or not sex plays a role in an individual bird’s movement after a territorial playback.

lmmtit <- lmer(dist_m ~ sex + (1|sex/id), data = tit)
summary(lmmtit)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: dist_m ~ sex + (1 | sex/id)
##    Data: tit
## 
## REML criterion at convergence: 220.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.3793 -0.5409 -0.0890  0.4035  1.6695 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  id:sex   (Intercept) 258.264  16.071  
##  sex      (Intercept)   4.187   2.046  
##  Residual              99.407   9.970  
## Number of obs: 28, groups:  id:sex, 14; sex, 2
## 
## Fixed effects:
##             Estimate Std. Error     df t value Pr(>|t|)    
## (Intercept)   28.838      6.533 14.604   4.414 0.000534 ***
## sex2          -2.629      9.909 14.236  -0.265 0.794568    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##      (Intr)
## sex2 -0.659
## optimizer (nloptwrap) convergence code: 0 (OK)
## unable to evaluate scaled gradient
##  Hessian is numerically singular: parameters are not uniquely determined

Using a linear mixed model to see if there’s an interaction between the sex of the bird and distance of detection.

anova(lm(dist_m ~ sex/id, data = tit))

## Analysis of Variance Table
## 
## Response: dist_m
##           Df Sum Sq Mean Sq F value Pr(>F)
## sex        1   47.4   47.40  0.1309 0.7207
## sex:id     2   90.6   45.31  0.1251 0.8830
## Residuals 24 8692.3  362.18

\(~\) \(~\)